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ABSTRACT 

Global three dimensional magnetohydrodynamic (MHD) simulations of turbulent accretion disks 
are presented which start from fully equilibrium initial conditions in which the magnetic forces are 
accounted for and the induction equation is satisfied. The local linear theory of the magnetorotational 
instability (MRI) is used as a predictor of the growth of magnetic field perturbations in the global 
simulations. The linear growth estimates and global simulations diverge when non-linear motions - 
perhaps triggered by the onset of turbulence - upset the velocity perturbations used to excite the 
MRI. The saturated state is found to be independent of the initially excited MRI mode, showing that 
once the disk has expelled the initially net flux field and settled into quasi-periodic oscillations in the 
toroidal magnetic flux, the dynamo cycle regulates the global saturation stress level. Furthermore, 
time-averaged measures of converged turbulence, such as the ratio of magnetic energies, are found to 
be in agreement with previous works. In particular, the globally averaged stress normalized to the gas 
pressure, < ap > = 0.034, with notably higher values achieved for simulations with higher azimuthal 
resolution. Supplementary tests are performed using different numerical algorithms and resolutions. 
Convergence with resolution during the initial linear MRI growth phase is found for 23 — 35 cells per 
scalcheight (in the vertical direction). 

Subject headings: accretion, accretion disks - magnetohydrodynamics - instabilities - turbulence 



1. INTRODUCTION 

Accretion disks are ubiquitous in astrophysics and 
play an essential part in the formation of stars and 
galaxies. For accretion through a disk to be effec- 
tive, angular momentum must be transported radially 
outwards, allowing material to move radially inwards. 
One means of achieving this is through viscous torques 
(|Lvnden-Bell fc Pringle!ll974D , and considerable progress 
has been ma de using the phenomenolog ical a-viscosity 
introduced bv lShakura fc Sunvaevl (I1973T ) which assumes 
that viscosity is provided by turbulent stresses which 
scale with the gas pressure. However, despite its suc- 
cesses, the a-viscosity model provides little physical in- 
sight into the mechanism(s) responsible for the turbulent 
stress. 

Even prior to the work of IShakura fc Sunvaevl (fl97l . 
instabilities in magnetize d rotat ing pla smas had been 
discov ered by IVelikhovl (|1959f ) and iChandrasekharl 
(fl960l). Yet, it was not until the seminal work of 
IBalbus fc Hawlevl (11991L 119981) that the so-called mag- 
netorotational instability (MRI) received widespread at- 
tention as the agent responsible for the onset of accretion 
disk turbulence. Linear stability analysis has shown that 
the MRI will amplify a seed magnetic field indefinitely 
until c onfronted by the strong-field limit or the diffusion 
scale (IBalbus fc Hawlevl 119921 iTerquem fc Papaloizoul 
119961: iPapaloizou fc Terqueml 119971 ). Non-linear stabil- 
ity analysis finds that growth of the magnetic field by 
the linear phase of the MRI is likely to be truncated 
by saturation re sulting from secondary, or pa rasitic, 
instabilities (e.g. iGoodman fc Xul fl99i iPessahl [2010h . 
That saturation of the magnetic field occurs was clearly 
demonstrated by even the very first shearing box sim- 
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ulations ([Brandenburg et ail 119951 : lHawlev et al.l 119951 : 
IStone et al.lll996fl . 

Contemplating the next steps in magnetized disks 
studies is aided by summarising what we have al- 
ready learned. For example, as mentioned above, 
it is clear that the magnetic field reaches saturation 
and that the resulting Maxwell stress dominates the 
angular momentum transport. In numerical simula- 
tions this necessitates high resolution to ensure that 
the fastest growin g MRI modes are sufficiently well re 
solved (see, e.g., ISanq et al.l 120041: IFromang fc Nelson 

120 



2551 lNoble"eTaLll2010t iFlock et allfeoilt lHawlev et al 



20111) . Related to this point is the importance of strat- 



ification, which introduces a characteristic length scale, 
removing the problem of non-convergence with simula- 
tion resolution encountered in unstratified simulations 
(IFromang fc Papaloizoul 120071: ILesur fc Longaretti|[2007 



Simon et al.l 120091: iGuan et al.l l2009t iDavis et alIl201C ; 
Sorathia et al.l I2012T T Stratification could also play a 



role i n the dynamo process which sets the saturation 
stress (lBrandenburg|[2005t lVishniad[200l iShi et alJl20Tol : 
iGressell l2010T l . However, the shearing box approxi- 
mation used in a large num ber of numerical studies 
to-date has limi t ations (e.g. iRegev fc Umurhanl 120081 : 
iBodo et al.l 120081 120111) . including the use of shearing- 
periodic boundary conditions in the radial direction, 
and/or periodic boundary conditions in the vertical di- 
rection. There boundary conditions artificially trap mag- 
netic flux, assisting the maintenance of the turbulent 
dynamo and obscuring the dependence of the saturated 
state on resolution. This is supported by a compari- 
son of peri odic and open boundary c onditions in global 
models by IFromang fc N elson ( 2006j) where the former 
were found to assist the dynamo by preventing magnetic 
flux from being expelled from the domain. In this re- 
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gard global models have the advantage of removing the 
vmphysical influence of the shearing box boundary con- 
ditions, albeit at a much larger computational expense. 

Other motivations for global models are the results 
from stability analyses of non-axisymmetric disturbances 
in magnetized accretion disks, where the most robust 
MRI modes are localized and the most robust buoyant 
(Parker) modes are global (iTerouem fc Papaloizoull 19961 : 
iPapaloizou fe Terqueml |1997|) . Therefore, large radial 
extents are required to accommodate the more global 
modes, and in this regard there is a limit to the ra- 
dial periodicity adopted in most shearing box simula- 
tions. These factors point to the need for high resolution, 
global, stratified disk simulations to further unravel the 
complexities of magnctorotational turbulence. 

Of the global simulation studies that have been 
performed a large number of the findings from lo- 
cal models have been maintained or have persisted; 
the ratio of the Maxwell and Reynolds stress is 
~ 3, and variations in toroidal magnetic fie ld with 
time are suggestive of a dynamo cycle (IHawlevI 



lHawlev fe Krolild 120021: iFromang fe Nelsonf 



iLvra et al 



Flock et al.l [20101 



120081: ISorathia et al.1 120101 
l201lL 120121 JO'Neill etaD 



2006 



2012 



2011 



Hawley et alj|2011t iBeckwith et al.H20lH iMignone et al 



2012t iMcKinnev et al.l l2012t iRomanova et all I2012D . 

However, a large number of these simulations do not 
start from fully equilibrium initial conditions where the 
magnetic field is accounted for in the force balance 
and in the induction equation. Both local and global 
models started with poloidal fields which do not sat- 
isfy the induction equation show rapid disruption and 
re-arran geme nt of the disk (e.g. iMiller fc Stoni 120001 : 
IHawlevI 120001 lHawlev et al]l201lh . This introduces a 
transient phase where channel flows are fueled by rapid 
shearing of the poloidal field lines. As such, extended 
run times are required to ensure that transients have 
subsided. To our knowledge, no previous global simu- 
lations of the MRI in stratified disks have used a fully 
equilibrium initial disk (i.e. satisfying force balance and 
the induction equation). 

We aim to explore the influence of magnetic fields on 
an accretion disk with global simulations. In this first 
paper we present equilibrium initial disk models with ar- 
bitrary radial density and temperature profiles. We then 
investigate the saturation (both locally and globally) of 
the growth of magnetic field perturbations. To this end 
we excite the MRI in global simulations using linear MRI 
calculations as a guide, and recover growth of magnetic 
field perturbations in agreement with estimates. In so 
doing we show that non-linear gas motions saturate the 
initial growth of the magnetic field and that at later times 
the turbulent state retains no knowledge of the initially 
excited MRI modc(s). 

The plan of this paper is as follows: in § [5] we de- 
scribe the equilibrium initial conditions and details of 
the numerical calculations and in §[3] we perform a linear 
perturbation analysis for the non-axisymmetric MRI. We 
present a suite of global magnetized disk simulations in 
§ |4j which explore the effect of different MRI mode exci- 
tation and numerical algorithms. In §[5] we compare our 
results to previous work, and then close with conclusions 
in §H 



2. THE MODEL 

2.1. Simulation code 

For our global disk simulations, we use a 3D spherical 
(r, 0, <f>) coordinate system with a do main which closely 
encap sulates the initial disk (e.g. IFromang fc Nelsonl 
120061 ) , and we solve the time-depe ndent equations of idea l 
MHD using the PLUTO code (|Mignone et all 120071 ). 
We note that throughout this work we describe our re- 
sults in terms of both spherical (r, 9, </>) and/or cylindri- 
cal (i?, <f), z) coordinates, with R = r sin and z — r cos 6. 
The relevant equations for mass, momentum, energy con- 
servation, and magnetic field induction are: 



| + V.[pv]= 0) 
^ + V ■ [pvv - BB + (P + i|B| 2 )I] = - P V$, 



at 



(i) 

(2) 



^ 771 

_+V-[(25 + P)v-(v.B)B] 



-pv ■ V$ - pA 

(3) 

— =Vx(vxB)(.4) 

Here E = pe + ^p|v| 2 + i|B| 2 , is the total gas energy 
density, e is the internal energy per unit mass, v is the 
gas velocity, p is the mass density, and P is the pressure. 
We use the ideal gas equation of state, pe = P/("y — 1), 
where the adiabatic index 7 = 5/3. The adopted scal- 
ings for density, velocity, temperature, and length are, 
respectively, 

/Oscalc = 1-67 x 10" 7 gm s" 1 , 
Vo = c, 

T scalc = pmc 2 /k B = 6.5 x 10 12 K, 
Z sca i e = l-48x 10 13 cm, 

where c is the speed of light, and the value of Z S caic cor- 
responds to the gravitational radius of a 10 s M black 
hole. 

The gravitational potential, $ of a central point 
mass (ignoring self-gravity of the disk), 3> is modelled 
using the pseudo-Newton ian potential introduced by 
iPaczvnskv fc Wiital (fl980l) : 

* = — V (5) 

r — 2 

Note that we take the gravitational radius (in scaled 
units), r g = 1. The Schwarzschild radius, r s = 2 for 
a spherical black hole and the innermost stable circular 
orbit (ISCO) lies at r = 6. The A term on the RHS of 
Eq ((3]) is an ad-hoc cooling term used to keep the scale- 
height of the disk approximately constant throughout the 
simulations; without any explicit cooling in conjunction 
with an adiabatic equation of state, dissipation of mag- 
netic and kinetic energy leads to an increase in gas pres- 
sure and, consequently, the disk scaleheight over time. 
The form of A is particularly simple, 



A 



P 



T(R,z)-T (R) 



(7 - 1) 2irR/ V4> 



(6) 



where Tq(R) and T(R, z) are the position dependent ini- 
tial and current temperature, respectively, is the ro- 
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tational velocity, and R is the cylindrical radius. This 
cooling function drives the temperature distribution in 
the disk back towards the initial one over a timescale of 
an orbital period and i s similar in its purpo s e to the cool- 
ing fu nctio ns used by IShafee et al.l (|2008D . INoble et al.l 
(|2010ft . and lO'Neill et al.l (l201lD~ Note that we only ap- 
ply cooling within \z\ < 2H, where H is the scaleheight of 
the disk, allowing heating via dissipation to occur freely 
in the corona. Our choice of an orbital period for the 
cooling timescale is somewhat arbitrary but is chosen as 
it represents a characteristic timescale for the disk. 

The PLUTO code was configure d to use the five- 
wave HLLD Riemann solver of Mi voshi fc Kusanol 
(12005ft. piece-wise parabo lic reconstruction (PPM - 
iColella fc Woodward! I1984D . and second-order Rungc- 
Kutta time-stepping. In order to maintain the V • B = 
constraint for the magnetic field we us e the upwind Con- 
strain ed Transport (UCT) scheme of IGardiner k, Stone! 
(|2008ft . Such a configuration has been shown to be ef- 
fective in recovering the linear growth rates of the ax- 



isymmetric MRI bv lFlock et~aH (|2010ft . In §g3]we test 
a number of different numerical setups: order of recon- 
struction, slope limitcrs, and simulation resolution. How- 
ever, in all of the other global simulations presented in 
§E]we use reconst ruction on characteristic variables (e.g. 
iRider et all 12007ft . A Courant-Friedrichs-Lewy (CFL) 
value of 0.35 was used for all simulations. 

The grid used for the global simulations is uniform in 
the r and (f> directions and extends from r = 4 — 34 and 
<f> = — 7r/2. In the 9 direction we use a graded mesh 
which places slightly more than half of the cells within 
\z\ < 2H of the disk mid-plane with a uniform A9, where 
H is the thermal disk scaleheight, and the remainder of 
the cells on a stretched mesh between 2H < \z\ < 5H. 
For example, for simulation gbl-mlO the 256 cells in 
the 8 direction arc distributed so that 140 cells are uni- 
formly spaced between 2H < \z\ < 5H. The respective 
grid resolutions and number of cells per scaleheight for 
the three global simulations are noted in Table [2J The 
grid cell aspect ratio at the mid-plane of the disk and 
at a radius of r = 18.5 (i.e. the disk midpoint) are 
rA9 : Ar : rAcj) = 1 : 1.4 : 8.6 and 1 : 1.6 : 2.5 for 
models gbl-mlO and gbl-ml0+, respectively. The r and 
9 boundary conditions depend on whether the cell adja- 
cent to the boundary contains > 1% disk material - which 
we determine using a tracer variable. If this constraint 
is satisfied we use outflow boundary conditions on all 
hydrodynamic variables except v$ which is determined 
from a zero-shear boundary condition (i.e. dfl/dr = 0) 
and the normal velocity, for which we enforce zero inflow. 
If the condition on disk material at the boundary is not 
satisfied we use outflow boundary conditions on hydro- 
dynamic variables with the limit that the values must lie 
between the floor values and the initial conditions for the 
background atmosphere - we find this choice to be useful 
in setting up a steady background inflow during the early 
stages of the simulation before material initially in the 
disk evolves to fill the domain. For the magnetic field 
we use zero gradient boundary conditions on the tangen- 
tial field components and allow the UCT algorithm to 
calculate the normal component so as to satisfy the di- 
vergence free constraint, with the exception that at the 
inner radial bounda ry we enforce a negat ive magnetic 
stress condition (e.g. iStone fc Pr inglc 200l|) . A periodic 



boundary condition is used in the <f> direction. Finally, we 
use floor density and pressure values which scale linearly 
with radius and have values at the outer radial boundary 
of 10~ 4 and 5 x 10~ 9 , respectively. 

2.2. Initial conditions 

Motivated by the fact that magnetorotationally turbu- 
lent disks are dominated by toroidal field, we start from 
an analytic equilibrium disk with a purely toroidal mag- 
netic field. The disk equilibrium is derived in axisym- 
metric cylindrical coordinates (i?, z); further details can 
be found in Appendix [X] along with alternative disk solu- 
tions which may be of use in future work. In the follow- 
ing we briefly summarise the equations for the isothermal 
in height, T = T(R), constant ratio of gas-to-magnctic 
pressure, (3 = 2P/\B\ 2 = 2P/B 2 , net magnetic flux disk 
adopted for the simulations presented in this paper. The 
choice of temperature and magnetic field lead to a den- 
sity distribution, in scaled units, 



p(R, z) = p(R, 0) exp 



-{$(i?,z)-$(i?,0)} P 



T{R) 



1 + B 



(7) 

where the pressure, P = pT. For the radial profiles 
p( R, 0) and T(R) we u s e sim ple functions inspired by 
the iShakura fc Sunvaevl (|1973[) disk model, except with 
an additional truncation of the density profile at a spec- 
ified outer radius: 



p(R,0) = p f(R, Ro,R 

out ) 



(8) 
(9) 



where po sets the density scale, Rq and R ou t are the 
radius of the inner and outer disk edge, respectively, 
f{R, Ro, Rout) is a tapering function and is described in 
Appendix [A] and e and x set the slope of the density 
and temperature profiles, respectively. In all of the global 
simulations R = 7, R out = 30, p = 10, T = 4.5 x 10~ 4 , 
€ = —33/20, and x = —9/10 (consistent with the radial 
scaling in the gas pressure an d Thomson-scattering opac - 
ity dominated region from IShakura fc Sunvaevl 11973ft . 
producing disks with an aspect ratio, H/R = 0.05. The 
rotational velocity of the disk is close to Kcplerian, with a 
minor modification due to the gas and magnetic pressure 
gradients, 



v%R, z) = vl(R, 0) + {$(R, z) - 0)}- — , (10) 
where, 

,d$(R,0) 2T 



v^(R,0) = R- dR 
1 + /A / RT dp{R,0) 



P J \p(R,0) dR 



(11) 



One advantage using such an equilibrium disk is that 
one begins with a disk that is close to the expected scale 
height and density. An isothermal disk, for example, has 
a scale height that is proportional to i? 3 / 2 . 

Finally, the region outside of the disk is set to be an ini- 
tially stationary, spherically symmetric, hydrostatic at- 
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mosphcrc with a temperature and density given by, 



PatmO") =Patm('"rcf) 



$(r) 
$(r re f) 



(12) 
(13) 



where p atm = 4 x 10 _5 po and r re f is a reference radius 
which we take to be i? ma x, the radius of peak disk density 
(see Appendix IX)) . The transition between the disk and 
background atmosphere occurs where their total pres- 
sures balance. 
As an example, model gbl-mlO corresponds to a disk 



with a peak density of 1.67 x 10 ' 
temperature of 2.9 x 10 9 K. 



gm s and a peak 



2.3. Diagnostics 

Turbulence is by its very nature chaotic. Therefore, 
averaged quantities are particularly useful diagnostics. 
In this section we describe how we calculate averages, 
and define the variables used to analyse the simulations. 

To compute shell-averaged values (denoted by curly 
brackets) of a variable q at a radius r we average in the 
8 and d> directions via, 



{<?} = 



/ qr 2 sin 8d6d(j) 
J r 2 sin 8d8d<j) ' 



(14) 



Similarly, we calculate a horizontally averaged value (de- 
noted by square brackets) as, 



[?] 



/ qr sin 6drd(j) 
J r sin 9drd<j) 



(15) 



To attain a volume-averaged value (denoted by angled 
brackets) we integrate over the radial profile of shell- 
averaged values and normalize by the radial extent, 



< q >- 



IU}dr 
Jdr 



(16) 



Time averages receive an overbar, such that a volume 
and time averaged quantity would read < q >. (Note 
that density-weighted averages are computed, but only 
for hydrodynamical variables.) For the analysis pre- 
sented in § U] we restrict the integration over r and 8 
to the range 10 < r < 30 and in tt/2 — 8 2 h/r < 
8 < tt/2 + 6 2H /r, where 9 2 h/r = taaT 1 {2H/R). 
We define this region as the "disk body" and limit 
the integration over this re gion to allow comparison 
against recent_ global (e.g. iFromang fc Nelson! 120061: 



Beckwith et all 120111: ISorathia et all 120101: iFlock et alj 
20111 12012b iHawley et al.1 120111 : ISorathia et all 1201. 



and large locaQ si mulations (e.g. lGuan fc Gammidl20Tl 
ISimon et al.l 120121) . 

In order to keep a track of the fluctuations in the scale- 
height of the disk during the simulation - which results 
from the interplay between adiabatic heating and our 
cooling function - a density-weighted average disk scale- 
height is computed, where we take H/R = c s /v§ (where 



1 Guan & Gammic (2011) find that properties of turbulence 
within the disk body in stratified d isks are quantitatively similar to 
those of unstratified disks (sec also Hawlcv ct al. 1995; Stone ct al. 
1996). 



c s is the sound speed), then perform a density- weighted 
shell-average followed by a radial averaging to acquire a 
volume averaged value, < H/R >. 

For accretion to occur, angular momentum must be 
transported radially outwards by turbulent stresses, and 
a major focus of numerical simulations is quantifying the 
stress. To this end, we define the perturbed flow velocity 
aa3 Svi = Vi — J V[r sin9d(f>/ J r sinOdcf) with i = R, </>, and 
compute the R — 4> component of the combined Reynolds 
and Maxwell stress, 



WRtj, = pdvnSv^, - Bb.B^,, (17) 
which is normalized by the gas pressure to acquire, 

{W m } 



{a P } 



{P} 



(18) 



Furthermore, we calculate the R — <f> component of the 
Maxwell stress normalized by the magnetic pressure, 



{«m} 



{\B\ 2 } ■ 



(19) 



To examine the operation of dynamo activity in the 
disk we compute the toroidal magnetic flux, which is de- 
fined as, 



B^r sin 9drd9. 



(20) 



The ability of the simulations to resolve the fastest 
growing MRI mod es is quantified in the sam e fashion as 
iNoble et al.l ([20T0l) and lHawlev etaLl (|20T1 . The wave- 
length of the fastest growing MRI modes with respect to 
the grid resolution in the z and <f> directions are, respec- 
tively, 



Amri- 
Az 



16 2tt\vaz\i" sind 
v<f,Az 



and, 



Amri-</> _ 27t|va0| 



RAcj) 



Acj) 



(21) 



(22) 



where «Az and tja^ are the vertical and azimuthal 
Alfvcn speeds, respectively, Af? and Acj) are the cell 
sizes in the 8 and <p directions, respectively, and Az = 
\J (r sin 8A8) 2 + (Ar cos 8) 2 is the corresponding cell size 
in the z direction. We define a single valued measure 
of resolvability as the fraction of cells in the disk bod y 
(|jz| < 2H) that have Q > 8 (e.g. ISorathia et alll2012jl . 



SC(Qj > 8) 
EC 



(23) 



where i = z,<f> and C represents a cell. The principal 
aim of calculating N z and is to quantify how well re- 
solved the turbulent state is in a simulation, and conse- 
quently whether global simulations are approaching the 
region of convergence fou nd from shearing box simula- 
tions (jHawlev et al.ll2011| ). 

2 Using an azimuthally averaged velocity when calculating the 
perturbed velo city removes the i nfluence of strong vertical and ra- 
dial gradients jFl'ock et al.ll201U ). 
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TABLE 1 

Parameters used for the linear MRI growth calculations 



Model 


P 


m 


fez 


fcR 


(^) 


Svo/Ca 


lin-mlO 


20 


10 


5 


2.5 


0.03 


0.1 


lin-ml0-/3100 


100 


10 


5 


2.5 


0.006 


0.1 


lin-ml0-/3300 


300 


10 


5 


2.5 


0.002 


0.1 


lin-ml0-/3300s 


300 


10 


5 


2.5 


0.002 


0.3 


lin-m40 


20 


40 


80 


40 


0.45 


0.1 



2.4. Fourier analysis 

To allow a direct comparison between the growth of 
MRI modes estimated from a linear perturbation analysis 
(§ [3]) and the results of global simulations (§ we anal- 
yse the growth of magnetic field perturbations in Fourier 
space. The procedure we follow is to remap the disk 
body (which we define in § 12. 3[) to a cylindrical mesh with 
uniform cell spacing in all directions, and a sufficiently 
fine resolution to ensure that the smallest cells from the 
spherical simulation grid are sampled. We then perform 
a 3D Fourier Transform of the data on the cylindrical 
grid. A detailed description of the cylindrical Fourier 
transform can be found in Appendix [5] 

3. EXCITING THE MRI 

Given that our global simulations commence with an 
equilibrium disk the MRI requires a seed perturbation to 
excite the growth of the magnetic field and development 
of turbulence. For this purpose we have chosen to excite 
a specific Fourier mode of the MRI using poloidal velocity 
perturbations. In the following we present perturbation 
calculations for the local, linear, non-axisymmetric MRI, 
the results of which are used in §H]to elucidate the evolu- 
tion of magnetic field perturbations in global numerical 
simulations. 

3.1. Linear MRI growth models 

Studies of the linear, non-axisymmetric MRI in weakly 
magneti zed disks have been examined by a number of 
authors dBalbus fe rlawlevlll992t iTerauem fc Papaloizou 
119961 iPapaloizou fc Terqueml 119971 ). iBalbus fc Hawlev 
(119921 ) 's local study showed that even if the seed mag- 



netic field is purely toroidal then the instability is still 
present, albeit with growth rates roughly an order of 
magnit ude lower than those fo und for initially poloidal 
fields (|Balbus fe Hawlevl Il99l . This result' was sup- 
ported by growth timescales approaching an orbital pe- 
rio d (for certain parameters) in m ore-global calculations 
bv lTerquem fc Papaloizou (1996|) where radial gradients 
were preserved. Furthermore, these authors found that 
in the k z /k R 1 limit - the primary domain of the MRI - 
instabilities become increasingly localized with time. On 
the other hand, in the k R /k z <^ 1 limit the Parker insta- 
bility dominates. In fact, even in the presence of dissipa- 
tion, MRI modes continue to become increasingly local- 
ized over tim e due to the time dependence of the radial 
wavenumber (IPapaloizou fc Terquemlll997D . Common to 
these studies is the finding that the non-axisymmetric 
MRI acts as a mechanism for the transient amplification 
of seed magnetic/velocity field perturbations by many 
orders of magnitude over tens of orbits. One question 
is, how well does this immense field amplification carry 
through to global, fully non-linear simulations? To an- 



swer this one needs an estimate of the linear growth. In 
this regard our analysis of the non-axisymmetric MRI in 
this paper is complementary to studies of the axisymmet- 
ric MRI in previous sim ulations (e.g. IHawlev fc Balbusl 
[i^9ltlFlocket^[2TjiTl . 

To construct our prediction for the global simulations 
we ut ilize the linear MRI model of IBalbus fc Hawlevl 
(|1992f) . (The perturbation analysis used to quantify the 
linear MRI growth is performed in cylindrical coordi- 
nates (R,(f>,z), whereas the global models presented in 
§ |H a re performed in sp herical coordinates (r,6, </>).) In 
brief, [Balbus fc Hawlevl perform a linear stability analy- 
sis of a local patch of a disk using the shearin g-sheet ap- 
proximation dGoldreich fc Lvnden-Bell 119651 ) where the 
perturbations are assumed to have a spatial dependence 
exp[i(k R R + m<fi + k z z)]. The equations for the evolution 
of the magnetic field perturbations form a pair of coupled 
second-order ordinary differential equation^. We let N 
be the Brunt- Vaisala frequency, which for the equilib- 
rium disk described in § 12.21 is. 



*r2 2 1 

N = - — 

5T 



1 



1 + /3J \RJ (r-2) 4 ' 
and define an independent time variable, 

dfl 



k 



■ k R {t)R = k R {t = 
k R + ^1 + K- 



0)R 



d\nR 



(24) 

(25) 
(26) 



Replacing the angular velocity with that due to a 
Paczynski-Wiita potential in the thin disk limit (i.e. 
H/R < 1), i} 2 = 1/R 2 {R - 2), the equations describ- 
ing linear perturbations ar^l, 



d 2 SB 7 , 



1k 7i 



dT 2 



1)1 



Rk 2 (3R - 2) 
4 / R-2 x 
~ \3R-2 



2r 2 

— 2-(i?-2)-#-2 



SB* 



(k-v A ) 2 

n 2 

. 4k z T 



k 2 -kl 



R 



k 2 
- 2 



k 2 m 2 \3R-2 



dSB T . 

dT 

dSB, 



dT 



(27) 



d 2 SB R 

dT 2 

2 R-2 
R 2 k 2 3R-2 

m 



k 2 
ft R. 

k 2 



:{t 2 -R 2 V 



R 



3R-2 
R + 2 



2 , dSB, 
-Rk. 



' dT 

dSB R 



JL ( R ~ 

~ \ZR 



(k- va) 

n 2 



-5B T - 



R-2 
k^N^ 

Rk 2 n 2 



dT 
tSB 7 , 



(28) 



3 Note that there i s a ty pographical error in equation (2.19) 
of IBalb us & Hawlev ( 1992]) where the final term should read 
&BzN 2 (kl -k 2 )/k 2 . 

4 The angular velocity resulting from our disk model (cf Eqs ilOt 
and Ijllll ) actually includes a small offset to Keplerian rotation. 
However, we find that this makes little difference to the pertur- 
bation calculations, and the subsequent comparison against global 
simulations in § [4] Therefore, for the sake of simplicity, we adopt 
a purely Keplerian rotation profile for the local calculations. 
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where SB Z and SB R are the vertical and radial magnetic 
field perturbations. 

The time dependence of k R in Eq (|25|) is a consequence 
of the radial wavenumber be ing sheared. Th e refore , 
within the framework of the IBalbus fc Hawlevl (|1992[ ) 
analysis the radial wavenumber can grow indefinitely so 
that radial disturbances can evolve to arbitrarily small 
spatial extent. Clearly, when we come to making a com- 
parison against our global simulations, this will not be 
the case due to finite numerical resolution. 

The magnetic field perturbations are related through 
the divergence- free constraint, 



k ■ SB = k R SB T , 



R 



SBa + k z 5B z = 0. 



(29) 



The unperturbed magnetic field topology only enters 
through k ■ va- For our initially purely toroidal magnetic 
field one finds, 



n 



2(R - 2)m 2 T 



(30) 



To initiate the MRI we use the R and z components of 
the linearized induction equation, 



d5B R 
dt 



i(k • B)6v R , 



and, 



dSB z 
dt 



= i(k • B)Sv z , 



(31) 



(32) 



where Sv R and 6v z are the poloidal velocity perturba- 
tions (with the imaginary part of Sv corresponding to 
the real part of dSB/dt). For the perturbations in the 
z-components in both the linear MRI and global calcu- 
lations we use a waveform, 



8v z = 8vq cos(k R R + m<j) + k z z), 



(33) 



which, on substitution into Eq (|32[) , and with the con- 
version between real and imaginary parts accounted for 
by a phase shift in the trigonometric term, leads to, 

dSB z 2B R-2 . 

— — = dv Q — -sm(k R R + m(f) + k z z), (34) 

dt 6R—2 

where Svq is the amplitude of the initial velocity pertur- 
bations. An equivalent treatment to Eq (|34f is used for 
the perturbations in the i?-componcnts with the differ- 
ence that we make use of the incompressibility condition, 
k ■ Sv = 0, and set, 

k 

Sv R = Svq- — cos(k R R + m<j> + k z z), (35) 
k R 

The remaining parameters used in the calculations are 
summarised in Table [T] 

Our first calculation, model lin-mlO, uses a f3 = 20 
magnetic field and wavenumbers for the excited MRI 
mode of m = 10, k z = 5, and k R = 2.5. These 
wavenumbers are chosen to ensure sufficient resolution 
in the global simulations and we leave a more detailed 
discussion to § |H The amplitude of the initial velocity 
perturbations, Svq, is set to 0.1 c s , where c s (= VT) is the 
sound speed. Since we intend to use these calculations as 
a guide for our global simulations, we use the equilibrium 
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Fig. 1. — The evolution of magnetic field perturbations from 
linear MRI growth calculations showing successful magnetic field 
amplification. From top to bottom: 8B^,SB Z , and SB^. The 
parameters used in these calculations are provided in Table [T] 



disk model described in § !2.2l to choose the input density 
and temperature. Calculations are performed at a cylin- 
drical radius, R = 20, and at the disk mid-plane where 
iV 2 = (see Eq (25)). From Eq © the disk tempera- 
ture, T = 1.75 x 10~ 4 , and the density, p = 0.46. The 
initial components of SB are set to zero, so too is the ini- 
tial azimuthal velocity perturbation, Svs - the poloidal 
velocity perturbations seed the instability through the 
dSB/dt terms. To integrate Eqs ([27]) and (28} we use 
an adap tive stepsize, 4th -order, explicit Runge-Kutta 
method (Press et al.l 119861) . As Fig. [T] shows, the mag- 
netic field perturbations grow extremely rapidly over the 
first few with noticeable oscillations, where pp rb 

is the radially dependent orbital period of the disk at 
cylindrical radius j. The upper panel of Fig. [2] shows 
the effective j3 for the MRI mode - the time required for 
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Iin-m10 
Iin-m10-pi00 




c ' 1 1 1 1 1 1 1 1 1 

2 4 6 8 10 12 14 16 18 20 



Time (P^,°) 

Fig. 2. — (Top): The plasma- f} calculated from the linear MRI 
growth calculations - see Table [TJ for the list of models and pa- 
rameter values pertaining to each calculation. The vertical line 
indicates when model lin-mlO reaches /3 = 1. (Bottom): The time- 
dependent radial wavenumber. The red and blue lines corresponds 
to models lin-mlO and lin-m40, respectively. The horizontal lines 
indicate the fcp values at the Nyquist limit for the global simula- 
tions gbl-mlO and gbl-mlO-lllr (see Tables \3\ and [2] and §|4j. 

the magnetic field to grow to (3 = 1 is only a few or- 
bital periods for model lin-mlO. Evaluating the approxi- 
mate growth rate, w of the magnetic energy, (as the 
gas pressure remains constant) via f5~ x = (3q exp(wt), 
we find an average growth rate over the first six orbits, 
uj = 0.14 Q. Applying the same approach to 5-Br we 
find uj = 0.090. This is cons istent with the findings of 
iTerquem fc Papaloizoul (|1996l ) but is roughly an order of 
magnitude larger than values of a few percent of the or- 
bital frequency quoted in gene ral for the development o f 
the non-axisymmctric MRI bv lBalbus fc H awlev ( 19921 ). 
Keeping all parameters fixed and then varying the initial 
magnetic field strength, one sees from models lin-mlO- 
/3100 and lin-ml0-/3300 the trend that the growth rate 
of <5B decreases with increasing initial (3. In model lin- 
ml0-/3300s the size of the initial velocity perturbations 
is increased to 5vq = 0.3 c s with the result that over the 
very first few orbits the growth of <5£Ts becomes very sim- 
ilar to that of a stronger initial field strength excited by 
smaller velocity perturbations. For a higher wavenum- 
ber perturbation the rate of initial growth increases, as 
evidenced by model lin-m40 (see Figs. [JJ and [2]). Evalu- 
ating the approximate growth rate of the magnetic field 
energy and (5£>r, for model lin-m40 gives, U = 0.68 51 and 
0.25f2, respectively From these results one may predict 



0.2 




-0.1 - 



-0.2 - 



20 40 60 80 100 

Time (P^ b ) 

Fig. 3. — Evolution of <5-Br (upper), 6B Z (middle), and o m 
(lower) as a function of time in model lin-mlO over a longer time 
duration than shown in Fig. [l] 

that the development of £B in simulations will depend 
on the initial field strength and/or the wavenumber of 
the excited mode(s). In § I4.3I we examine if this result 
h olds true in global si mulati ons. 

IBalbus fc Hawlevl (Il993) discuss the parameter 
(1c-va) 2 /0 2 and attribute to it an important role 
in the ability of the MRI to successfully amplify the 
seed field. They find that for (k-v A ) 2 /^ 2 > 2.9 the 
instability is stabilized and magnetic field oscillations 
are damped. For the models shown in Fig. [TJ this 
parameter is much less than unity. From Eq (|30[) for 
(k-VA) 2 /^ 2 one can see that to increase the value of 
this variable one can either decrease (3 - which increases 
the tension along field lines - or employ higher azimuthal 
wavenumbers, m. The latter has the side-effect of 



For comparison, the maximum growth rate for the axisymmet- 



ric MRI is 0.75H IIBalbus fe HawlevlUMlD . 
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increasing the growth rate of fep and causing tight 
wave crest wrapping, both of which lead to a more 
rapid stabilization of the radial disturbances. However, 
we find that irrespective of the value of (k-VA) 2 /^ 2 , 
which is 0.03 for lin-mlO (Table [T]), the perturbation 
in the magnetic field ultimately decays. This is shown 
in Fig [3] (upper and middle panels) where the lin-mlO 
calculation is plotted for a longer time duration. Despite 
continuing growth in SB Z , there is decay in <5-Br, which 
is a consequence of the increase of k^(t) combined with 
the divergence-free constraint (Eq (|25]l). The ratio of 
the Maxwell stress to magnetic pressure (a m ) predicted 
from the linear MRI growth calculations is shown in the 
lower panel of Fig. [3j where we define, 




a m = 2SB *^ + B *\ B * = Bl + SBl + 5Bl + SBl 

(36) 

Clearly, considering that a m , or to be more exact 

< aM > (its global analogue - see Eq (fl9]) ), is 
a commonly used diagnostic in numerical simulations 
(jHawlev et al.|[l995l) . Under the action of the linear MRI 
alone < o?m > would never reach a steady value. This 
ulti mate decay of linear MRI distu rbances is consistent 
with lTerquem fc F^a paloizoul (|1996lV s finding of transient 
instability growth in a number of numerical tests, in 
which k z > fcpt initially. Shear causes &r to grow but 
once fcp{ > fc z growth halts. 

The ultimate decay of linear MRI modes has sig- 
nificance for global models because the maintenance 
of dynamo action requires all components of the field 
to be sufficiently strong. This highlights the need 
for an additional mechanism, other than the linear 
MRI growt h, to replenish SB^ (e.g. parasitic in- 
stabilities - iGqqdman fc Xul Il994|: Pa rker instability - 
iTout fc Pringld 119921 IVishniad 120091: dynamo action 
in the steady-state turbulence iBrandenburg et al.lll995l 
lHawlev et aflU996h . 

The linear MRI growth calculations act as a check on 
our global simulations, principally to examine whether 
our setup recovers the growth rates of the linear MRI 
accurately. However, there is a limit to the time interval 
when we can confidently make a comparison between the 
linear grow th models and global sim ulations. Firstly, the 
analysis of iBalbus &: Hawlevl (|1992l ) adopts the Boussi- 
nesq approximation which becomes invalid when the az- 
imuthal magnetic field becomes super-thermal. The up- 
per panel of Fig. [2] shows that this limit is reached in 
approximately 2.3 orbits for lin-mlO and 0.1 orbits for 
lin-m40. Secondly, fcpt(i) can grow indefinitely in the 
linear growth models, yet this is not the case for our 
global simulations which are restricted by finite numeri- 
cal resolution. Taking the Nyquist limit to be 2 grid cells, 
and considering, for example, the resolution of model 
gbl-mlO, the maximum resolvable radial wavenumber is 
fcjt_Nyquist = 86. This limit is reached after ~ 16.5 and 
2.3 orbits for models lin-mlO and lin-40, respectively 
(lower panel of Fig. [2]). Therefore, choosing to excite 
a higher wavenumber MRI mode limits the time interval 
where comparisons can be made against linear pertur- 
bation theory, and this is one reason why we choose to 
excite a lower wavenumber mode (m = 10) in the global 
simulations. 




2 4 6 8 10 12 14 16 18 20 22 24 26 28 



Time (P° 3 ' °) 

Fig. 4. — The evolution of < ap > (upper) and < Oiyi > (lower) 
in the global models as a function of time in units of the orbital 
period at R = 30. See Table [3] for a list of the models and time 
averaged values. 

4. GLOBAL MODELS 

In this section we describe the results of global simu- 
lations using the initial conditions and simulation setup 
described in § [3] The global simulations are listed along 
with grid dimensions, number of cells per scaleheight, 
and approximate MRI growth rates in [2j Time and vol- 
ume averaged variables quantifying the steady-state tur- 
bulence are given in Table [3] In models gbl-mlO and 
gbl-ml0+ we excite a specific Fourier mode using a plane 
wave, which takes the form of Eq (|33|) . as described in 
§ [3] These models use the same wavenumbers as model 
lin-mlO so as to allow a direct comparison of magnetic 
field growth. The third model, gbl-rand, uses random 
pressure and poloidal velocity perturbations to initially 
seed the disk disturbance. All of the global models start 
with a purely toroidal magnetic field with j3 = 20. Mod- 
els gbl-ml0+ and gbl-rand are computed on grids with 
lower poloidal resolution (roughly 2/3 that of model gbl- 
mlO), but with a factor of three better azimuthal resolu- 
tion. In the following section we present some properties 
of our model disks and demonstrate that higher <j) resolu- 
tion to be a crucial ingredient in producing a sustained, 
high valued turbulent stress, < ap >. 

4.1. Model evolution 

We begin with a description of the evolution of models 
gbl-mlO and gbl-ml0+ (Tables [5] and [3J. In this model 
we adopt an azimuthal wavenumber which varies with 
cylindrical radius, m — m(R). We give the azimuthal 
wavenumber a radial dependence of m{R) = m cr it(-R)/6, 
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TABLE 2 

Global simulations and the corresponding linear growth rate. 



Model 


Resolution 


Reconstruction 


Limitcr 


n r /H 


n e /H 


n<f,/H 


k-'approx 




(n r x n g x nj,) 






(10 < r < 30) 


(1*1 < 2H ) 


(n(R = 20)) 


lin-mlO 














0.09 ± 0.01 


gbl-mlO 


768 x 256 x 128 


Parabolic 


Char 


18-77 


35 


4 


0.15 ± 0.03 


gbl-mlO+ 


512 x 170 x 320 


Parabolic 


Char 


12-51 


27 


12.5 


0.11 ± 0.02 


gbl-rand 


512 x 170 x 320 


Parabolic 


Char 


12-51 


27 


12.5 


0.15 ± 0.04 


gbl-mlO-lin 


768 x 256 x 128 


Linear 


Char 


18-77 


35 


4 


0.16 ± 0.03 


gbl-mlO-cw 


768 x 256 x 128 


Parabolic 


CW84 


18-77 


35 


1 


0.13 ± 0.02 


gbl-mlO-cs 


768 x 256 x 128 


Parabolic 


CS08 


18-77 


35 


1 


0.14 ± 0.02 


gbl-mlO-hr 


896 x 300 x 150 


Parabolic 


Char 


21-90 


11 


5 


0.16 ± 0.04 


gbl-mlO-lr 


512 x 170 x 96 


Parabolic 


Char 


12-51 


23 


3 


0.13 ± 0.03 


gbl-mlO-llr 


342 x 112 x 64 


Parabolic 


Char 


8-34 


15 


2 


0.09 ± 0.04 


gbl-mlO-lllr 


192 x 64 x 32 


Parabolic 


Char 


5-19 


9 


1 


0.05 ± 0.03 



TABLE 3 

Time averaged quantities from the global simulations.. 



Model 


m 




A'R 


c 


5tav a 


N z 


iV,, 


< B l> 


<B|> 
<B'i> 


<B|> 

< B l> 


< ap > 


<«M> 


< /S > b 


gbl-mlO 


10 


5 


2.5 


0.1 


10-28 


0.32 


0.33 


0.071 


0.25 


0.018 


0.017 


0.31 


32 


gbl-ml0+ 


10 


5 


2.5 


0.001 


12-26 


0.45 


0.75 


0.12 


0.29 


0.036 


0.034 


0.41 


17 


gbl-rand 








0.001 


16-26 


0.45 


0.75 


0.13 


0.30 


0.037 


0.034 


0.41 


17 



Note. — a Time interval over which averaging was performed, 
b Time averaged plasma- /3 in the disk body. 




10 12 14 16 18 20 22 24 26 28 
Time (P§[f) 

Fig. 5. — Rcsolvability of the MRI in the cf> direction (upper) and 
z direction (lower) in models gbl-mlO, gb l-m l0+, a nd g bl-rand. 
For a definition of the resolvability see Eq l|23| l and § I2.3I 

where the criticaQ azimuthal wavenumber for the linear, 



rr 
x 




Fig. 6. — Density-weighted and volume averaged scalehcight of 
the disk, < H/ R > as a function of time in models gbl-mlO (solid 
line) and gbl-mlO-l- (dashed line). 



non-axisymmetric MRI, 



merit 



(37) 



At R = 20 in model gbl-mlO, m cr it = 60; adopting m = 
m cr it / 6 ensures that the corresp onding fc z and &r are wel l 
resolved by the numerical grid. iBalbus &: Hawlevl ()1992l ) 
noted that the fastest growing non-axisymmetric modes 



occur for fc z 



'/R, and we also use this relation to 



calculate fc z . Given that our grid resolution is coarser in 
r than it is in 9, we set = k z /2. 

The initial poloidal velocity perturbations seed the 
growth of magnetic field perturbations via the MRI and 
after roughly 1 — 2 P^ h turbulent motions become ap- 



6 Defined as the value of m for which distur bances grow most 
rapidly, which follows from equation (2.30) of IBalbus fc Hawlevl 



11992ft. The radial dep endence of Ec I (37} stems from the 
IPaczvnskv Wiital (flOSOf l potential (Eq (P). 
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parent in the disk body. As the poloidal magnetic field 
becomes established throughout the disk the resulting 
Maxwell stresses disrupt the disk equilibrium. The evo- 
lution of models gbl-mlO and gbl-ml0+ is largely similar 
during the first few orbits of the simulations. Examin- 
ing the normalized stress, < ap >, shows that there is 
an initial transient phase which peaks after a simulation 
time of roughly 4 P^ h (Fig. 0). Following this, < ap > 
gradually decreases until a steady-state is reached after 
roughly 12 P^ h and the time-averaged stress for the 
remainder of the simulation, < ap > = 0.017 for gbl- 
mlO. The time-averaged ratio of the Maxwell stress 
to the magnetic energy, < an > — 0.31, which is be- 
low the values of r oughly 0.4 quoted by, for example, 
lHawlev et al.l (|2011l ) for well resolved turbulence. To in- 
vestigate the dependence of these values on the azimuthal 
resolution of the simulation we have also run model gbl- 
ml0+, which has 12.5 cells/if in the <f> direction (and a 
lower resolution in the poloidal direction - see Tabled]). 
The higher <f> resolution clearly influences the turbulent 
stresses in the simulation and for model gbl-mlO we find 
< «m > = 0.41 and < ap > = 0.034, in agreement with 
high resolution shearing-box simulations. The resolvabil- 
ity of the fastest growing MRI modes (see Eq (|23[) and 
Fig. also clearly show that a higher azimuthal resolu- 
tion helps to maintain (or even strengthen) the poloidal 
magnetic field - models gbl-mlO and gbl-ml0+ initially 
show similar values of N z but largely different values of 
N^, and combined with the evidence mentioned above 
is evident that azimuthal resolution is very important 
for maintain ing a healthy turbulent state (see also the 
discussion in 



lHawlev et al 



Fromang fc Nelsonll2006l; iFlock et al.ll2011t 
|2011[ ). In ?[5l we compare further quantita- 



tive measures of the steady-state turbulence to previous 
works. 

The poloidal magnetic field develops in flux tubes with 
small spatial scale, which dissipate magnetic energy via 
reconnection, heating the disk. In Fig. [5] we show the 
density- weighted and volume-averaged scaleheight of the 
disk, < H/R > as a function of time. In model gbl- 
mlO the scaleheight of the disk increases initially un- 
til t ~ 8 P^ h , after which it steadily declines. This 
shows that during the initial disk evolution, dissipation 
heats the disk more rapidly than the cooling function, A 
(see Eq ([6])) can drive the temperature back to its ini- 
tial value. In other words, the dissipative timescale is 
shorter than an orbital period. In contrast, for model 
gbl-ml0+, < H/R > remains roughly constant after the 
initial rise, which shows that the higher < ap > in this 
model (Fig.01) is causing more heating, and a marginally 
thicker disk. 

Fig. [7] shows a poloidal slice through the disk in model 
gbl-ml0+ at t = 14 P§q ■ During the turbulent steady 
state the disk is characterised by a dense, cold, sub- 
thermally magnetized core close to the mid-plane and 
a tenuous, hot, trans-to-super thermal magnetic field at 
z > 2H (the corona). Turbulent motions are clearly 
evident in the plot of in Fig. [7] with the dominant 
eddies appearing to have a larger size in the corona com- 

ared to the disk body. As noted by iFromang &; Nelson] 



20061) . such behaviour arises due to conservation of an- 
gular momentum in eddie motions - or wave action - as 
small scale eddies rise out of the dense disk mid-plane 



■ 
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Fig. 7. — Slices in the poloidal plane from model gbl-mlO-t- at 
t = 14 P^q° showing p (upper), T (middle), and (lower). 
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Fig. 8. — Disk body, corona, and global volume averages of the 
plasma-/? for model gbl-ml0+. 



Exciting MRI modes in equilibrium disks 



11 



1e-02 



1e-03 - 



m 1e-04 



1e-05 



1e-06 




1000 



Fig. 9. — Comparison of |<5_Br (»?i = 10,fc z = 5, kn = 2.5) | from 
models lin-mlO (local linear MRI) and models gbl-mlO, gbl-ml0+, 
and gbl-rand (fully non-linear global simulations). 



into the less dense coronal region. Our intended pur- 
pose for the explicit cooling function, A (see § I2.ll for 
details) becomes more apparent from the temperature 
plot - we aim to take a step beyond the purely isother- 
mal approximation and towards the obscrvationally sup- 
ported picture of a hot corona and cooler disk body. 
In Fig. [8] we show the volume-averaged plasma-/?. In 
the disk body, we find < (3 > = 17 and for the corona 
6. The coronal value is higher than values of 

Stond 



</3> 

close to one found in pr evious isothermal (j Miller 
[2000HFlock et al.ll20ll an d quasi-isothermal sim ulations 
(jFromang fc Nelson! I200l iBeckwith et~all 1201 1[ ). which 
may be attributable to the lack of any explicit cooling 
in the corona in our simulations. However, although 
the gas in the corona is heated by dissipation, it does 
not continually heat up through the simulation, and in 
fact remains quasi-steady through the latter half of the 
simulation. This contrasts with adiabatic shearing-box 
simulations with imposed periodic bo undary cond i tions , 
in which the ga s does heat up (e.g. iStone et"aT1 119961 : 
iSano et ahl 120041 ) and demonstrates that when coronal 
gas is allowed to freely expand, adiabatic cooling can, to 
some extent, balance heating via turbulent dissipation. 

4.2. Comparison with linear MRI growth estimates 

In Fig. |9] we compare the evolution of |<5Br| for the 
m = 10, fc z = 5, &r = 2.5 mode (measured in Fourier 
space - see § I2.4|) for model lin-mlO (which describes lin- 
ear MRI growth - see Tables [2] and [T]) and models gbl-mlO 
and gbl-ml0+ (global simulations which allow fully non- 
linear evolution - sec Table [3]) . Quantifying the initial 
growth by deriving approximate growth ratefl, w a pprox 
for the curves shown in Fig. [9] , we find that the linear 
MRI estimate is matched best by model gbl-ml0+, with 
gbl-mlO (and gbl-rand) producing higher growth rates. 
The higher amplitude perturbation for model gbl-mlO 
compared to gbl-TOl0+ originates from the larger am- 
plitude initial velocity perturbation of 0.1 c s (compared 
to 0.001 c s for gbl-ml0+ and gbl-rand - see Table [3]). 
The agreement between the global simulations and linear 
growth estimate (model lin-mlO) begins to falter after 
roughly 2 P^q and growth in |<5-Br| for the global simu- 

~ An approximate exponential growth rate is determined by fit- 
ting \SB R \ = (5-Br|o exp(w a p prox t). 
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Fig. 10. — Plasma-/3 from models lin-mlO and gbl-mlO (upper), 
and time evolution of some sample non-linear terms in model gbl- 
m.10 (lower). 

lations levels off. From the linear MRI calculations shown 
in Fig. [2] one may anticipate that the growth in \SB\ is 
halted by the magnetic pressure evolving to equipartition 
with the gas pressure - which is illustrated by the vertical 
dashed line in Fig. [5]- and would mean that one cannot 
rely on model lin-mlO as a predictor for model gbl-mlO. 
However, Fig. 1101 shows that j3 remains roughly constant 
in model gbl-mlO over the first few orbits. What then 
causes the local MRI estimates and the global simula- 
tions to diverge? In the lower panel of Fig. [TU] we plot 
the evolution of the following non-linear terms derived 
from the momentum equations, 



dv r 
1 dr 

dr 



b = 



,d = 



r 86' 
ve dv d 
r 89 



Small spikes in these non-linear terms occur after ~ 
0.1 orbit. However, after 1-2 orbits considerably larger 
fluctuations become apparent, particularly in b which 
also appears to have the highest amplitude oscillations 
of the plotted terms thereafter. Considering that the 
linear MRI is seeded by velocity perturbations through 
the induction equation (Eqs (f3Tj) and ([32]) ), the corre- 
lation in time between the non-linear velocity terms be- 
coming active and the growth in |<5.Br| departing from 
the linear MRI growth predictions is highly suggestive 
of non-linear motions causing saturation in the growth 
of a specific MRI mode. Furthermore, the small am- 
plitude kicks from these non- linear terms after 0.1 or- 
bits may explain the early divergence between the /3 
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1e-09 



1e-10 



Fig. 11. — Logarithmic false-colour image showing the time evo- 
lution of |5Ba(k)| in model gbl-mlO+. Values of 5Br in the disk 
body were Fourier transformed and results are shown at specific 
values of m and k z , and the full range of fe^. From top to bottom: 
m = 10 and k z = 5 (low wavenumber), m = 40 and fc z = 80 (mod- 
erate wavenumber), and m = 120 and k z = 150 (high wavenum- 
ber). 



values predicted from model lin-mlO and those found 
from gbl-mlO. In this sense the non-linear motions pro- 
vide saturation to the initial phase of local SH growth. 
Whether the non-linear motions are attributable to sec- 
ondary insta bilities feeding off the linear MRI growth 



onclary insta bilities feeding on tne linear ivltu growtn 
locally (e.g. iGoodman fc Xul H99l IPessah et al J [20071: 
iPessah fc Goodman! 120091 IPessah! 120101) . or are due to 
the onset of turbulence ( Latter et al.l 120091 ) propagating 
radially outwards through the disk is unclear and would 
require an analysis of the non-linear growth of the non- 
axisymmetric MRI, which we do not pursue here. 

In summary, comparisons between linear growth cal- 
culations and global simulations highlights a number of 
potential saturation mechanisms. Such as, growth of 
magnetic field perturbations beyond the weak field limit, 
and/or growth of the radial wavenumber beyond the fi- 
nite limit of the simulation resolution. However, for the 
simulations performed in this work, we find that satura- 
tion of growth in magnetic field perturbations correlates 
well with the onset of non-linear motions. 
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Fig. 12. — Same as Fig. [TT] except for model gbl-rand. 

4.3. Trigger dependence 

A major focus of magnetized accretion disk simulations 
is to study properties of the quasi-steady-state turbu- 
lence. A necessary test is whether the turbulent steady 
state depends on the MRI mode initially excited, and also 
whether prohibitive transient behaviour arises due to the 
choice of exciting a specific MRI mode. For this purpose 
we have computed model gbl-rand which uses the same 
initial disk as model gbl-ml0+ with the difference that 
the disk is perturbed with random perturbations in the 
both the poloidal velocity (amplitude 5vo = 0.001 c s ) 
and gas pressure (10% amplitude) which excite a range 
of MRI modes. Simulation resolution and time-averaged 
measures of the turbulent state are listed in Tables [2] and 
[31 respectively. 

The evolution of model gbl-rand is very similar to that 
of model gbl-ml0+; the initial perturbations excite the 
MRI and lead to growth of (5-Br. Both models show sim- 
ilar growth in the m = 10, k z , /cr. = 2.5 mode (Fig. O 
which one would expect given that this mode is excited 
with the same amplitude perturbation. After 3 
(~ 0.5 P§Q h ) values of 8Bn become almost identical be- 
tween the models irrespective of the differing initial per- 
turbations. We illustrate this in Figs.[TT1and[T2lin which 
we show the evolution of SBr in Fourier space for a range 
of /cr values. The different panels in the figures corre- 
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spond to low, moderate, and high wavenumber values for 
m and k z (relative to the size of the disk and the Nyquist 
limit). As mentioned above, 5 Br values are very similar 
between the two models at t > 3 P£g h . Furthermore, 
even though we excite a specific low wavenumber mode 
in model gbl-mlO, a wide range of modes rapidly emerge. 
We attribute this behaviour to wave-wave coupling and 
the onset of a turbulent cascade. 

Exciting larger wavenumbers should provide a larger 
initial MRI growth rate (see § [3]), but how does this af- 
fect the evolution of magnetic field perturbations in the 
global simulations? In particular, does the wavenumber 
of the excited MRI mode affect the globally-averaged 
saturation stress? In model gbl-rand a white noise 
spectrum of perturbations has been excited. Therefore 
higher wavenumber modes can contribute to the ini- 
tial growth phase in < ap >. There is an indication 
of this from Fig. ITS"]) where the growth of \8B-r\ at a 
range of wavenumbers means that the Maxwell stress, 
and consequently < ap > will also grow across a range 
of wavenumbers. Fig. [4] shows that < ap > does grow 
faster for gbl-rand compared to gbl-mlO+ (which have 
identical grid resolution) , supporting the notion that the 
growth in the globally averaged stress due to an ensemble 
of unstable modes is higher than for a single wavenumber 
mode. 

All three models start with a toroidal magnetic field 
with a net flux, and during the early evolution of the 
disk, the combination of magnetic buoyancy and accre- 
tion expels magnetic flux from the disk body such that 
by the time the turbulent steady state is reached the net 
toroidal magnetic flux of the disk, is close to zero. 
Subsequently, vE^ oscillates about the zero-point with a 
period of roughly 5 orbits (upper panel of Fig. [13]) consis- 
tent with previo us global simulations and suggestive of 
a dynamo cycle dFromang fc Nelson! [20061 : 1 O'Neill etaD 
120111: iBeckwith et all 120111 ). All three models demon- 
strate this behaviour. However, minor differences in the 
toroidal magnetic flux, vE^, are visible between models 
gbl-ral0+ and gbl-rand (Fig. H3]). The different models 
are slightly out of phase, which is not surprising given 
the differences in the transient evolution at early simu- 
lation times (Fig. HJ . Interestingly, model gbl-rand does 
not overshoot when expelling the initial net toroidal flux 
and thus settles into dynamo oscillations at a slightly 
earlier time, which may explain why the transient phase 
in < ap > takes a longer time to fade in this model. 

In conclusion, once the disk reaches a turbulent steady- 
state the disk retains no knowledge of the MRI mode ini- 
tially excited. This is supported by the almost identical 
time- averaged properties of the disk noted in Table [3] for 
models perturbed by a single low wavenumber mode or 
an ensemble of modes. 

4.4. Algorithm and resolution dependence 

In this section we examine the ability of different nu- 
merical algorithms to recover the growth of magnetic 
field perturbations resulting from the non-axisymmetric 
MRI. Comparisons between numerical simulations and 
analytical esti mates for the axisymmetr ic M RI have been 
presen ted by lHawlev fc Balbua (|1991l) and iFlock et all 
(|2010D . Considering that MHD turbulence m accre- 
tion disks produces a predominantly toroidal field it is 
important to examine how well numerical algorithms 
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Fig. 13. — Toroidal magnetic flux at (j> = 0, ^^((p = 0) as a 
func tion of time for models gbl-mlO, gbl-ml0+, and gbl-rand. See 
Eqs (l2"0l for the definition of 

can recover the growth of <5B as a result of the non- 
axisymmetric MRI. The setups used are listed in Table [5J 
The different combinations are intended to test different 
orders of reconstruction, parabolic limiters, and grid res- 
olution^]- Reconstruction refers to the order of accuracy 
used to interpolate cell interface values (which are then 
used in the Riemann solver to calculate fluxes of con- 
served variables) . Parabolic limiters are used to preserve 
monotonicity and prevent extrema from being intro- 
duced by the reconstruction step. We examine the orig- 
inal li miter for PPM proposed by iColella fc Woodward! 
(119841) . the extremum p reserving limiters presented by 
IColella fc Sekoral (|2008l ). and limiters base d on recon- 
struc tion via characteristic variables (e.g. iRider et"all 
120071 ). The aforementioned slope limiters are respec- 
tively denoted "CW84", "CS08",and "Char" in Table [1 
The last parameter we vary is the grid resolution, as this 
places a constraint on the maximum resolvable wavenum- 
ber, and for this purpose we compute models gbl-mlO- 
hr, gbl-rolO, gbl-rnlO-lr, gbl-mlO-Ur, and gbl-mlO-lllr 
(which have decreasing resolution). Note that with the 
exception of models gbl-ml0+ and gbl-rand, all models 
have the same cell aspect ratio and the same ratio of cells 
in the disk body to cells in the corona as gbl-mlO. 

As in § 14.21 we calculate approximate growth rates of 
the magnetic field perturbation, \SBr (k)| for the m = 10, 
fc z = 5, kn = 2.5 mode via a Fourier analysis of the initial 
simulation evolution. The results are shown in Table [2] 
which we summarise as follows: 

• The best agreement with the linear MRI growth 
rate comes from model gbl-ml0+, showing that 
azimuthal resolution is important for properly cap- 
turing MRI growth. 

• Within errors the choice of limiter does not make a 
considerable difference to the resulting growth rate. 

• linear reconstruction produces a comparable 
growth rate to parabolic reconstruction. 

8 Our aim is to examine the ability of d ifferent algorithms to 
capture waveforms. Balsara & Meyer (2010) found that Riemann 
solvers that do not resolve the Alfven wave are more likely to lead 
to turbulence dying out as a result of a higher level of numerical 
dissipation. Therefore, we have not endeavour ed to test different 
Riema nn solvers and we use the HLLD solver of[Miyoshi & Kusanol 
(2005) in all global simulations. 
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• There is a consistent trend that growth rates in- 
crease with increasing resolution (see models gbl- 
mlO-lllr, gbl-mlO-lh, gbl-mlO-lr, gbl-mlO, and 
gbl-m 10-hr). With 4 cclls/_ff in the <f> direction 
the growth rates are converged (within errors) for 
roughly 23-35 cells/i? in the z direction. This is 
a slightly lower threshold than the > 40 ce\\s/H 
(in the vertical direction) found to achieve conver- 
gence in the time-averaged turb ulent stress in strat - 
ified shearing-box simulations (jDavis et al.ll2010T ). 
Comparison with model gbl-ml0+ (which has 12.5 
cells/ff in the (j> direction) suggests that conver- 
gence in global disks may be achieved at lower res- 
olutions when the cell aspect ratio is closer to unity. 

5. DISCUSSION 

With a growing number of studies using strati- 
fied shearing boxes with high re s olution and /or a 
large spatial extent fShi et all 120101: iDavis et all 120101: 
iGuan fc Gammid 120111: ISimon et all 120111 I2012D and 
higher resolution global models (IFromang fe Nelson! 
200(1 ISorathia et all 201(1 l20li iFlock et alJl201lL I2012t 



Beckwith et al. 1 12011b IHawleT et al .1 1201 lU Mi gnone et al.l 
20121 ) . quantifying the steady state turbulence and mak- 



ing direct comparisons between simulations provides a 
check of both consistency and convergence. 

One of the most popular measures of the steady state 
is < ap >. In this regard, models gbl-TOl0+ and gbl- 
rand produce valu es of ~ 0.034 wh i ch is higher than re- 
cently reported bv iBeckwith et al~l (|2011[ ) and, as noted 
by these authors, higher than previous global models 
and a num ber of high resolut ion shearing-box simu- 
lations (see lHawlev et al.l 1201 lL and references there- 
in). We attribute the larger < ap > in our mod- 
els to a higher azim uthal resolution than used by 
([Beckwith et aljfeoill ) , but also note the possible indica- 
tion that higher < ap >'s may be more readily achievable 
in global simulations. Our average < q;m > ~ 0.41 (for 
models gbl-ml0+ and gbl-rand) is in good agreement 
with the ~ 0.36 — 0.4 achieved by t he highest resolu- 
tion shearing-box s imula tions to-date (jDavis etalll20Tol: 
ISimon et all 1201 U I2012I ). Considering that our models 



have a lower number of cells/i? than the aforementioned 
shearing-box models, there may also be an indication 
that convergence may be achieved at lower grid resolu- 
tions than in localized models, potentially due to averag- 
ing over a larger volume, and capturing lower wavenum- 
ber eddies. 

Comparing models gbl-mlO and gbl-ml0+, strong evi- 
dence points to the grid cell aspect ratio and, in particu- 
lar, the resolution in the <j) direction as an important pa- 
rameter in ac hieving a high "< ap > and < qm > (see the 
discussion in Fromang fc Nelsonll2006l Flock et al.ll201lL 
lHawlev et all 120111 and ISorathia et al.ll2012h . A possi- 
ble explanation for this is that the dynamo cycle - which 
helps to sustain the turbulent state and involves the MRI 
as a driving agent - can operate more effectively at higher 
frequencies when the cell aspect ratio is closer to unity. 

Related to < aM > is the tilt ang l e, 8tiit, where 
sin_20 t ii t =< aM > (jGuan et al.l 120091: IBeckwith et all 
120 111) . It has been argued bv ISorathia et al.l (|2012ft 
that this parameter provides a better measure of con- 
vergence than < ap >, at least in the case of unstrat- 
ified turbulence for which the question of convergence 



in the absence of exp l icit d issipation was raised by 
iFromang fc Papaloizoul pOOl . We find 6 t iit ~ 9° for 
model gbl-mlO which is consiste nt with previou s find- 
ings for stratified global disks (IBeckwith et al.1 120111 
lHawlev et all 120 111 IFlock et al.H2012f ). However, models 
gbl-TOl0+ and gbl-rand have Otut ~ 12° which is compa- 
rable to values of ~ 11° — 13° for sheari ng box simulations 
(both unstratified and stratified, e.g.. IGuan et al.l 120091 : 
ISimon et al.ll2012|) and also for recent stratified global 
disks calculat ions performed with an orbital advection 
algorithm by iMignone et al.l (|2012| ). These results are 
encouraging as they show that global simulations are 
reaching sufficient grid resolution to reproduce shearing- 
box results. 

The ratio of directional magnetic energy also pro- 
vides insight into convergence and correspondence be- 
tween simulations 



Wc find, 



< B 2 

~ R 



>/< 



> 



> 



0.13, 
0.036 



< B 2 > I < Bl > ~ 0.30 and < B 2 > / < B 2 
for models gbl-m 10+ an d gbl-rand. These values are 
higher than obtained by lHawlev et al.l (1201 1| ) for their 
global disk simulations, and in some cases only slightly 
lower than values found from high r esolution sh e aring- 



box simulations ( Shi et all 120101: 
Simon et al.l 120111 : IGuan fc Gammiel 



20llr 



Davis et al.1 [2oTob 
20111: ISimon etall 



Interestingly, model gbl-mlO produces a sustained 
stress, albeit with a lower value than model gbl- 
ml0+, but with only 4 cells/if in the (^-direction. 
IFlock et al.l (|2011l) found that at least 8 cells/if were re- 
quired to produce a sust ained turbulent stress (see also 
IFromang fc Ne lson 2006). However, these authors used 
linear reconstruction, whereas we have used parabolic 
reconstruction which may permit a sustained stress at a 
slightly lower resolution. 

Finally, we note that we do not see any prominent ev- 
idence of recurring transient phenomena due to linear 
growth revi vals in the mean m agnetic fields, as recently 
reported bv IFlock et all poll . This may be due to dif- 
ferences in the nume rical setup between our models and 
those of IFlock et all or perhaps this phenomena occurs 
at later times that we have not reached with the simula- 
tion runtimes of our models. 

6. CONCLUSIONS 

We have performed global 3D MHD simulations of tur- 
bulent accretion disks which start from fully equilibrium 
MHD initial conditions. The local linear theory of the 
MRI is used as a predictor of the growth of magnetic field 
perturbations in the global simulations. Additional tests 
have also been performed to compare results obtained 
from global simulations using a number of different nu- 
merical algorithms and resolutions to the linear growth 
estimates. Our main findings are: 

i) The growth of magnetic field perturbations in the 
global simulations shows good agreement with the 
linear MRI growth estimates during approximately 
the first orbit of the disk. Subsequently, the over- 
whelming influence of non- linear motions, which may 
be due either to the onset of turbulence or to sec- 
ondary instabilities, saturates the local growth. 

ii) The saturated state is found to be independent of 
the initially excited MRI mode, showing that once 
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the disk has expelled the initial net flux field and 
settled into quasi-periodic oscillations in the toroidal 
magnetic flux, the dynamo cycle regulates the global 
saturation stress level. Furthermore, time-averaged 
measures of quasi-steady turbulence are found to be 
in agreement with previous work. In particular, the 
time averaged stress, < ap > ~ 0.034. 

iii) We find <B'l> I <B\> ~ 0.13 for global strat- 
ified simulations with 12.5 ce\\s/H in the cf> direc- 
tion, which is in good agreement with value found 
from high resolution, stratified, shearing-box simu- 
lations. Higher <fi resolution in the simulation (at 
least > 4 cells/if) is required to maintain stronger 
radial and vertical magnetic field, and consequently 
a larger < ap >. 

iv) From the numerical algorithms that we tested, the 
choice of reconstruction order or limiter does not 
significantly alter the resulting linear MRI growth 



rate. Convergence with resolution (for the linear 
MRI growth tests) is found for resolutions of roughly 
23 — 35 cells per scalehcight (in the vertical direc- 
tion). However, above all, a higher azimuthal reso- 
lution contributes to a much better agreement with 
linear growth estimates, supporting the push for low 
cell aspect ratio (close to one) in global accretion 
disk simulations. 
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APPENDIX 
A. EQUILIBRIUM DISK SOLUTIONS 

The initial conditions for our simulations are an equilibrium thin disk with a purely toriodal magnetic field. In the 
following we present some analytic solutions which are of use for numerical simulations of accretion disks and for disks 
in other environments, such as starburst galaxies (e.g. ICooper et al.ll2008l ). These solutions incorporate more or less 
arbitrary radial profiles of density and temperature and a toroidal magnetic field. The latter involves either a constant 
ratio of gas-to-magnctic pressure, /3, radially dependent /?, constant B^, and variants with net/zero toroidal magnetic 
flux. 

In axisymmetric cylindrical coordinates (R,(j>,z), in steady-state, and with = v z = -Br = B % = 0, the induction 
equation is identically satisfied and we are left with the two momentum equations: 

dP 9$ P v% 1 d(R 2 B 2 ) 

dR +P dR ~R~ + 2IT- dR ~ U ' (Ai) 

op a$ idB 2 , 

^ + ^ + 2^°< (A2) 

where the pressure, P = pT, with T in scaled units. 

We derive a compatibility condition for the above equations by subtracting d/dR of Eq (|A2|) from d/dz of Eq (|A1[) . 
to obtain, 

dT\_dp_ 1 dp d(R 2 B 2 ) i dB 2 dTldp 1 dpdBl 
dz p dR 2R 2 p 2 dz dR P R dz dRpdz 2p 2 dR dz [ ' 

To solve for the disk equilibrium we take the approach of using Eq (|A2[) to derive an equation for p(R,z), Eq (|A3[) 
to acquire v 2 (R,z), and Eq (|A1[) to obtain an expression for v^(R, 0). The resulting equations require boundary 
conditions for the run of p and T at the disk midplane, which can be chosen arbitrarily. From here on we take the disk 
to be isothermal in height, T = T(R), and firstly consider a disk with a constant, /3 = 2P/\B\ 2 = IP/B 2 ^. Eq (|A"2j) 
then becomes, 

l°p = _?*(_L\J_ (A4) 

pdz dz \l + pj T{R) { ' 

which integrates to give an expression for the density in terms of its midplane value, 




Turning to the rotational velocity, the compatibility relation (Eq (|A3p ) reduces to, 

'l + P\ IdpdT 




P J p dz dR 



(A6) 



which upon integrating and using Eq (|A4I) leads to an expression for the azimuthal velocity in terms of its midplane 
value, 

v 2 (R, z) = v 2 (R, 0) + {*(JJ, z) <f>(R, 0)}- m (A7) 

The model is completed with a midplane rotational velocity, which is determined by substituting Eq (|A5I) into Eq (| Al|) . 
This gives 

2 , B „, „ mR,0) , 2T (l + P\( RT dp(R,0) dT\ 
v,(R, 0) = R^- + -+ {— J [ p{m dR +R-). (A8) 

The first term is the square of the Keplerian velocity; the remaining terms are proportional to the square of the sound 
speed so that Eq (|A8[) represents a minor departure from a Keplerian disk. 

A possible variation to the aforementioned disk would be to make P radially dependent, i.e., P = P(R) . For example, 
one may choose to make P(R) oc sin(fci?), where k is a radial wavenumber. In such a case Eq (|A5j) for p(R,z) is 
unchanged. However, the expression for the rotational velocity becomes, 

vl(R, z) = vl(R, 0) + {*(*, ,) - $(R, 0)} - ^) 1) (A9) 

where, following substitution in Eq (|A1|) . 

2^ R r^ „ ^(i?,0) , 2T , (1 + P\( RT dp(R,0) , n dT\ RT dp 
v,{R, 0) = R^^- + - + [— ) [ p{R Q) QR +RjR)-^dR- (A10) 
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Alternatively, one may desire a disk with a constant B^,(R, z) = -B^o, in which case the magnetic pressure does not 
influence the density profile, leading to, 

P[R, z) = P\R, 0) exp I I, (All) 



and a corresponding velocity profile of, 

R dT 

v%R, z) = vl(R, 0) + z) - 0)}- — + vi^R, z) - v%{R, 0), (A12) 

with, 

2 , . <9$(i?,0) a . , / i?T dp(R,0) dT\ fl s 

and where the Alfven speed, vai/, (R,z)^B^/^pjR~zj. 

As mentioned above, the radial profiles p(R, 0) and T(R) required to complete the disk model may be chosen 
arb itrarily, subject to boundar y constraints at the outer disk edge. As an example, we use simple functions inspired by 
the IShakura fe Sunvaevl (|1973[) disk model, modified by truncation of the density profile at a specified outer radius: 

p(R,0) = p f(R,R ,R out )(^-) , (A14) 




f(R, Ro, R out ) = hR^ + t^-t =^ - 1 , (A15) 



t{r)=t °{r\T' (A16) 

where po sets the density scale, i?o and i? ut are the radius of the inner and outer disk edge, respectively, and e and 
\ set the slope of the density and temperature profiles, respectively. The tapering function, f(R, Ro, Rout) is used to 
truncate the disk at an inner and outer radius. In practice this function is normalized to give p(i? m ax,0) = po, where 
the radius of peak density, i? ma x is given by the positive root of the quadratic resulting from taking d/dR of Eq (IA14I) . 
namely, 

V^W = ^(VR^t + VR^) + ^^{^R^+^Rvf - 4(1 - e- 1 )^^ (A17) 

where a = 1 — (2e) . Once i? max is known it is straightforward to renormalize the density profile. 

Fi nally, studies of turbulent dynamos in magnet i zed disks are often concer ned with the net flux of the magnetic field 
(e.g. iBrandenburg et alJll995t lHawlev et al.l ll996: Fro mang fc N elson 20QU). For the initially purely toroidal field we 
have adopted in this paper the net flux of the disk is given by = J J B^dRdz. Noting that in the above derivations 

we have used f3 to relate B^ to P, meaning B$ = ±y/2pT/ (3, i.e. we are free to choose the sign of B^. Therefore, if a 
net flux field is required then one may set the sign of B^ the same everywhere, whereas if one desires a zero-net flux 
field then, for example, one may choose to make B^ anti-symmetric about the disk midplane. 

B. FOURIER TRANSFORM IN CYLINDRICAL COORDINATES 

We wish to evaluate the Fourier transform F(k) of a function /(r) = f(R, <fi, z) expressed in terms of cylindrical 
polar coordinates (i?, <f>, z). The definition of the Fourier transform is 

F(k) = J exp(ik • r)/(r) d 3 x (Bl) 

Cylindrical coordinates in real and Fourier space are expressed via the following equations: 

X = R cos 4> fc x = /cr cos ip 

y = Rs'mcf) k y = fcji sin tp (B2) 
z — - z k<z, ~ — 

Hence, 

k ■ r = k^R cos(cj) — -0) + k z z (B3) 

and 

F(k) = / exp[i(fc R i? cos(0 - ip) + k z z)] f(R, </>, z) RdRd^dz (B4) 
Jv 

where V is the computational region, usually of the form: 

R < R < Ri < cj) < 2tt - z < z < z (B5) 
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We begin by constructing a Fourier series in the periodic azimuthal coordinate </>: 

oo 

f(R,4>,z)= fm(R,z)e- im " f ' 

m— — oo 

where the coefficients f m (R, z) arc given by: 

f m (R, z) = -L J 2 J f(R, <t>, z y m * d<p 



(B6) 



(B7) 



We now make the change of angular variable x — 4> ~ V'! the integration over % is still over the interval [0, 2n] since all 
of the angular functions within the integrand have period 2-7T. The Fourier transform can now be expressed as: 



F <M= E 



Jmip 



e ik **f m (R,z) 



2tt 



i(-mx+k R R cos x) 



dz 



RdR 



The angular integral can be expressed in terms of Bcsscl functions (J m (k R R)): 

r 2n 



i(-mx+k n Rcos x) 



d X = 2iri m J m (k R R) 



n 



Hence, 



00 j-Ri 

F(k) = 2tt ^ i m e lm ^ I 

m=-oo •* R « 



J m (k R R) 



e lk * z f m (R,z)dz 



RdR 



(B8) 

(B9) 
(BIO) 



Equation (|B10p defines the following procedure: 



im<j> 



1. Evaluate the angular coefficients: 

i r 2 * 

f m (R,z) = — J f(R,<f>,zy 

2. Then perform the integration in the z direction: 

e ik * z f m (R,z) dz 

-zo 

3. Finally, perform the (truncated) Hankel transform in the radial direction: 

rRi 

F m (k R ,k z ) = / RJ m (k R R)F m {R,k z )RdR 

Jro 

4. The Fourier transform of f(R, <fi, z) is 

oo 

F(k R ,iP, k z ) = 2Tr i m e im ^ F m (k R ,k z ) 



(Bll) 



(B12) 



(B13) 



(B14) 



m— — oo 



Since the input data for /(i?, ^, z) are on a grid, the azimuthal, vertical and radial wave numbers, m, k z and k R , 
are limited by the Nyquist limit. Let the number of intervals in each coordinate direction be (n R , n^, n z ) and the grid 
increments be (Ai?, A</>, Az) = [(i?i — Ro)/n R , 2n /n^, 2zo/n z ]. The grid coordinates are <f> v , z w where: 



R u = Rq + uAR u = 0, 1, n R — 1 
</>„ = u A0 u = 0, 1, - 1 

z u , = — zq + w Az w = 0, 1, n z — 1 



(B15) 



The expressions for the azimuthal ,f m {R, z) and vertical F m (R, k z ) parts of the Fourier transform can be approximated 
by discrete Fourier transforms as follows: 



f m {R u ,Z w ) ps — V" f(Ru,(j> v ,z w ) cxp[2irimv/7 



"0 n 
r u= 



0, 1,..., n /2 



(B16) 



F m (R u ,k z ) ps — exp[-ifc z z ] f m (Ru, z w ) exp[27rZu;/n z ] i = 0, 1,..., n z /2 fc z = — / 



to=0 
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The radial transform can be evaluated from 



-2 



F m (kR, k z ) = k R / sJ m (s)F m (s/k R ,k z ) ds (B17) 

Jku_Ro 

More accurate versions of equations (|B16[) may be evaluated using the approach given in iPrcss et al.l (|1986[) . 



